function [ZeroRates, CurveDates] = zbtprice(Bonds, Prices, Settle, ...
    OutputCompounding, varargin)
%ZBTPRICE Bootstrapped spot and forward curve from Bond Prices.
%
%   [ZeroRates, CurveDates] = zbtprice(Bonds, Prices, Settle)
%   [ZeroRates, CurveDates] = zbtprice(Bonds, Prices, Settle,
%         OutputCompounding)
%
%   Inputs:
%    Bonds - The portfolio of coupon bonds from which the zero curve will
%            be derived, specifically, an NxM matrix of bond parameters
%            where each row of the matrix corresponds to an individual bond
%            and each column corresponds to a particular parameter.
%
%            The required columns (parameters) for this matrix are:
%              Maturity - [Column 1] the maturity for each bond in the
%                         portfolio in serial date number form
%
%            CouponRate - [Column 2] the coupon rate for each bond in the
%                         portfolio in decimal form
%
%            Optional columns (parameters) are:
%                    Face - [Column 3] Nominal value of bond. Default is
%                           $100.
%
%                  Period - [Column 4] the number of coupon payments per
%                           year in integer form.
%
%                           Possible values are:
%                           0,1, 2 (default), 3, 4, 6, and 12
%
%                   Basis - [Column 5] values specifying the basis for each
%                           bond in the portfolio.
%
%                           Possible values:
%                           0 - actual/actual(default)
%                           1 - 30/360 (SIA)
%                           2 - actual/360
%                           3 - actual/365
%                           4 - 30/360 (PSA)
%                           5 - 30/360 (ISDA)
%                           6 - 30/360 European
%                           7 - actual/365 Japanese
%                           8 - actual/actual ISMA
%                           9 - actual/360 ISMA
%                          10 - actual/365 ISMA
%                          11 - 30/360 ISMA
%                          12 - actual/365 ISDA
%                          13 - bus/252
%
%            EndMonthRule - [Column 6] value specifying whether or not the
%                           "end of month rule" is in effect for each bond
%                           contained in the portfolio.
%                           1 - on (default)
%                           0 - off
%
%   Prices - [Nx1] column vector containing clean price values per $100
%            face for each bond contained in the portfolio represented by
%            the Bonds matrix.
%
%   Settle - [Scalar] value representing time zero in derivation of the
%            zero curve. Normally this is also the settlement date for the
%            bonds contained in the portfolio from which the zero curve
%            will be derived.
%
%   Optional Inputs:
%   OutputCompounding - scalar value representing the period by which the
%                       output zero rates will be compounded; the default
%                       value is semi-annual (i.e. "2") compounding.
%
%                       Possible values are:
%                       1, 2, 3, 4, 6, 12, -1.
%
%   Outputs:
%    ZeroRates - Nx1 vector containing the values for the implied zero
%                rates for each point along the investment horizon defined
%                by a maturity date.
%
%   CurveDates - Nx1 vector containing the maturity date for each zero rate
%                along the investment horizon (from time T = Settle to time
%                T = maturity of the longest dated bond in the source
%                portfolio).%
%
%   Notes: 1) In cases where the source portfolio of bonds contains more
%             than one bond with the same maturity date, the mean zero rate
%             is calculated for that maturity date.
%          2) Ensuring that the source portfolio contains a sufficient
%             number of bonds and that those bonds are evenly distributed
%             with respect to maturity date will significantly enhance the
%             performance of this function.

% Copyright 1995-2006 The MathWorks, Inc.
% $Revision: 1.11.2.13 $   $Date: 2009/09/23 14:01:20 $

%Check to ensure that the minimum number of arguments has been passed in
if (nargin < 3)
    error('Finance:zbtprice:tooFewInputs', ...
        'The Bonds matrix, Prices vector and settlement date must be entered.');
end

[NumBonds, NumCols] = size(Bonds);
if (NumCols < 2)
    error('Finance:zbtprice:tooFewColumns', ...
        'The Bonds matrix must contain at least maturity date and coupon rate.')
elseif (NumCols > 6)
    error('Finance:zbtprice:tooManyColumns', ...
        'The Bonds matrix may only contain 6 columns.')
end

if (nargin > 4 || ~isempty(varargin))
    warning('finance:zbtprice:noLongerSupported',...
        ['Additional inputs are now ignored.  Please type "help ',...
        'zbtprice" for updated input arguments.']);
end

% Set default columns: Face 100, Period 2, Basis 0, EOMRule 1.
Col = ones(NumBonds,1);
DefaultCols = [100*Col, 2*Col, 0*Col, 1*Col];
Bonds = [Bonds, DefaultCols(:, (NumCols-1):4)];

% Parse out bond parameters
Maturity = Bonds(:,1);
CouponRate  = Bonds(:,2);
Face     = Bonds(:,3);
Period   = Bonds(:,4);
Basis    = Bonds(:,5);
EndMonthRule  = Bonds(:,6);

% Determine the compounding frequency
if all(isisma(Basis))
    CompFreq = 1;
elseif all(~isisma(Basis))
    CompFreq = 2;
else
    error('Finance:zbtyield:bondsInconsistentBasis', ...
        'All Bonds must have either all ISMA or SIA bases.');
end

% Prices
Prices = Prices(:);
if length(Prices)~=NumBonds
    error('Finance:zbtprice:pricesDoNotMatchBonds', ...
        ['The number of Prices, %d, does not match the ',...
        'number of bonds %d'],length(Prices), NumBonds)
end

% Parse Settlement
Settle = finargdate(Settle);
if any(Settle ~= Settle(1))
    error('Finance:zbtprice:bondsMustSettleOnSameDay', ...
        'All Bonds must Settle on the same day.');
else
    Settle = Settle(1);
end

% Set OutputCompounding
if nargin<4 || isempty(OutputCompounding)
    OutputCompounding = 2;
end

% Create the cash flow amounts, dates, and time factors
[CFAmounts, CFDates, CFTimes] = ...
    cfamounts(CouponRate, Settle, Maturity, ...
    Period, Basis, EndMonthRule, [], [], [], [], Face);

NumInst = length(Prices);

CFAmounts = [-Prices              CFAmounts];
CFDates =   [Settle*ones(NumInst,1) CFDates];
CFTimes =   [zeros(NumInst,1)       CFTimes];

% Lay out the cash flows by instrument and by date
[CFSet, Dates, Times] = cfport(CFAmounts, CFDates, CFTimes);
NumTimes = length(Times);

% Find the maturity of each cash flow stream (max date of nonzero cash
% flow)
[Maturity, I] = max(CFDates .* (CFAmounts~=0) ,[],2);

% Remember the maturity cash flow and time factor
CFEnd = CFAmounts( (1:NumInst)' + NumInst*(I-1) );
TFEnd = CFTimes(   (1:NumInst)' + NumInst*(I-1) );

% Find the ordered unique maturities for the output rates
[EndDates, I, InstOrder] = unique(Maturity);
NumPoints = length(EndDates);

% Make sure the problem is not overdetermined:
% Collapse cash flow streams which have the same maturity to reduce the
% NumInst cash flows to NumPoints synthetic cash flows.
[JPoint, InstOrder] = ndgrid(1:NumPoints, InstOrder);
MatMap = (JPoint == InstOrder);

% Reduce the rows of CFSet by summing cash flows.
CFSet = MatMap*CFSet;
CFEnd = MatMap*CFEnd;

% Just pick up any examples of unique maturities
TFEnd = TFEnd(I);

[CurveDates, EndInd] = intersect(Dates, EndDates);
EndTimes = Times(EndInd);

% Cheat with interp to do the work for now
% It could be faster to do this by hand
RateMap = zeros(NumTimes, NumPoints);

% Extrapolate rates before the first maturity (EndDate)
RateMap(1:EndInd(1)-1,1) = 1;

% loop over the points to get their linear interpolation weights
% Only interpolate times after the first maturity at index EndInd(1)
for j=1:NumPoints,
    EndRates = zeros(NumPoints,1);
    EndRates(j) = 1;
    RateMap(EndInd(1):end,j) = interp1q(EndTimes, EndRates, Times(EndInd(1):end));
end

% Initial guess for the rates ignoring intermediate cash flows
EndRates = 2*( (-CFEnd./CFSet(:,1)).^(1./TFEnd) - 1 );
EndRates = max(EndRates, 0.04);

% Objective function
    function [Error, dEidRj] = zerobootsub(EndRates)
        
        Rates = RateMap * EndRates;
        Disc    = (1 + Rates./CompFreq).^(-Times);
        DelDisc = (1 + Rates./CompFreq).^(-Times-1) .* (-Times./CompFreq);
        DelDisc(1) = 0;
        
        Error  = CFSet*Disc;
        dEidRj = CFSet*(DelDisc(:,ones(1,size(RateMap,2))).*RateMap);
    end

fsolveOpt = optimset('Jacobian','on', 'Display','off', ...
    'TolFun',1e-12,  'TolX',1e-12, ...
    'MaxIter', 20);

[ZeroRates, ~, ExitFlag] = fsolve(@zerobootsub,EndRates,fsolveOpt);

if ( ExitFlag <= 0 )
    warning('Finance:zerobootcf:solutionConvergenceFailure', ...
        'Could not solve for the rate\n');
end

% Transform to a different OutputCompounding if requested
if OutputCompounding~=2
    
    if OutputCompounding == -1
        ZeroRates = CompFreq*log(1 + ZeroRates/CompFreq);
    else
        ZeroRates = OutputCompounding * ...
            ((1 + ZeroRates/CompFreq).^(CompFreq/OutputCompounding) - 1);
    end
end
end